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ABSTRACT 

Flexion-based weak gravitational lensing analysis is proving to be a useful adjunct 
to traditional shear-based techniques. As flexion arises from gradients across an im- 
age, analytic and numerical techniques are required to investigate flexion predictions 
for extended image/source pairs. Using the Schwarzschild lens model, we demonstrate 
that the ray-bundle method for gravitational lensing can be used to accurately recover 
second flexion, and is consistent with recovery of zero first flexion. Using lens plane 
to source plane bundle propagation, we find that second flexion can be recovered with 
an error no worse than 1% for bundle radii smaller than AO = 0.01#e and lens plane 
impact pararameters greater than Oe + A0, where Oe is the angular Einstein radius. 
Using source plane to lens plane bundle propagation, we demonstrate the existence of 
a preferred flexion zone. For images at radii closer to the lens than the inner boundary 
of this zone, indicative of the true strong lensing regime, the flexion formalism should 
be used with caution (errors greater than 5% for extended image/source pairs). We 
also define a shear zone boundary, beyond which image shapes are essentially indis- 
tinguishable from ellipses (1% error in ellipticity). While suggestive that a traditional 
weak lensing analysis is satisfactory beyond this boundary, a potentially detectable 
non-zero flexion signal remains. 
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1 INTRODUCTION 

Weak gravitational lensing provides one of the most direct 
probes of the matter distribution of the Universe as it is 
independent of both the dynamical state and the nature 
of the matt er. Building on the pion e ering attempts by 
IValdes et aD <|l983h and iTvson et all (| 19841 ) to measure 
coherent changes in the shapes of background galaxies due 
to a foreground lens population, weak lensing techniques 
have now come of age. In recent years, weak lensing has 
successfu lly been applied in the cases of gala x y-galaxy lens- 
ing (e.g. Brainerd, Blandford & Smail 1996; Hudson et al 



2008) , lensing by clusters ( e.g. I Smail et~al 
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The conventional mathematical basis for weak lensing 
analysis assumes that a shear field causes an additional el- 
lipticity to the shape of a background source, which can be 
calcu lated by measuring the moments of the images ([Kaiserl 
1995). However, this approach does not account for higher- 
order shape distortions that occur when there are strong 
tidal fields across the image. Recently, an additional lensing 
effect called flexion has been investigated as an extension 
to shear-based measurements ([Goldberg &; Nataraiani r2002: 
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iGoldberg fe Baconl 120051 : iBacon et all l2006l fl. Flexion has 
two components denoted first and second flexion, which are 
a shift of the image centroid and a representation of the 
"arciness" of the image respectively. 

One main issue associated with shear-based gravita- 
tional lensing is that galaxies are intrinsically elliptical 
in shape (ellipses at some inclination to the line-of-sight 
seen in projection), so it is necessary to disentangle lens- 
induced shear from the intrinsic shape. Resolved galaxies, 
however, are not intrisincally flexed - although systems un- 
dergoing a merger, or galaxies with substantial asymmet- 
ric sub-structure such as a large starforming region, may 
be mis-intepreted as flexion signals. It has been suggested 
that fle xion may provide a stronger constraint on dark 
matter (ILeonard et al l |2007l: IBacon. Amara fe Readl l2010l : 
lHawken Bridle! 1 20091: iLeonard Kindl2010h. ga laxy clus- 
ter mass models ([Leonard. King fc W ilkins 2009) and also 
on delensing gravitational wave signals ([Shapiro et al.ll2010l ) 
than shear on its own, notwithstanding the challenges in 
measuring flexion (e.g. Okura, Umetsu & Futamase 2007, 
2008; Goldberg & Leonard 2007; Massey et al. 2007; Irwin 
& Shmakova 2006; Schneider & Er 2008). 

In a previous paper we presented analytic flexion results 
for a range of popular mass density profiles: Schwarzschild 
lens, singular isothermal sphere (SIS), Navarro- Frenk- White 
(NFW) profile and Sersic-like profiles (Lasky & Fluke 2009; 
hereafter Paper I). Our analytic solutions present a flexion 
formalism where we consider a two-dimensional (2D) field 
in the lens plane, which allows for the treatement of ex- 
tended sources. In this paper we extend our previous work 
by considering the following key question: over what range 
of image/source sizes is the flexion approach valid? That 
is, if flexion appears as a gradient of shear across an im- 
age, which is implictly assumed to be zero for traditional 
weak- lensing analysis, how well can we recover flexion for 
extended sources? 

The approach we use is via the ray-bundle method intro- 
duced by Fluke, Webster & Mortlock (1999). Here, bundles 
of light rays with a known initial configuration are propa- 
gated through one or more lens planes, and the deflection 
of each light ray is determined using the gravitational lens 
equation, see equation (pQ) below. The initial and final shapes 
of the bundles can now be used to obtain numerical esti- 
mates for convergence, shear, magnification, first and second 
flexion as a function of source size. By starting with simple 
lens models, we use known analytic results to test the accu- 
racy of our approach, and hence assess its applicability to 
cases where analytic flexion results do not exist (e.g. asym- 
metric mass profiles, such as from iV-body, dark matter halo 
simulations) . 

The paper is set out as follows: in Section [21 we sum- 
marise the key results for the analytic flexion formalism. 
In Section [3l we describe how the ray-bundle method can 
be used to obtain flexion along individual lines-of-sight. We 
present the results of our numerical testing, and determine a 
range of image locations and bundle radii for which the flex- 



1 see also [lrwii^& Shmakoval (l2Q05l . l20Q6h : 

llrwin. Shmakova fc Ande^onl (|2007h who independently de- 
rived equivalent higher-order gravitational lensing effects which 
they call sextupole lensing. 



ion approximation is valid for a Schwarzschild lens model. 
We demonstrate the existence of a flexion zone, which de- 
fines a physical region where the flexion formalism is valid 
for extended sources. In Section [H we present calculations of 
the size of the flexion zone as a function of lens mass, lens 
and source redshift. We present our conclusions and identify 
future directions for this work in Section [5] 



2 THE FLEXION FORMALISM 

In this section, we summarise relevant results of the flex- 
ion formalism - for full details, see Paper I and references 
therein. 



2.1 The gravitational lens equation 

The thin-lens gravitational lens equation relates the location 
of a background source and its image (s) due to an interven- 
ing mass distribution: 



Ds c n - 



(i) 



For the simple case of a single lens at the origin of the lens 
plane, is the image location in the image plane, r\i is the 
source location in the source plane, and on is the deflection 
angle. Throughout this paper, indices k = 1,2 signify 
the two-dimensional components, with summation assumed 
over repeated indices, and tensors are implicity functions of 
the image position, unless otherwise indicated. The thin- 
lens approximation requires that the spatial extent of the 
lens distribution be much smaller than the angular diame- 
ter distances between the observer and lens, Dl, observer 
and source, Ds, and lens to source, Dls- A convenient scal- 
ing of the lens equation into angular coordinates uses the 
substitutions /3i = rji/Ds, 0% — &/Dl and on — ^j-oti, so 
we have the compact form: 



ft 



(2) 



We introduce a further scaling for the lens equation in sec- 
tion 3.4. 



2.2 The lens matrix 

As a mapping between the lens and source planes, equation 
J2} is often rewritten in terms of its Jacobian matrix: 



dOi 



ft — 71 _ 72 
-72 1 - ft + 71 



(3) 



where the matrix elements are identified with the conver- 
gence, ft, and two orthogonal components of shear, 71 and 
72. The magnification along a given line-of-sight, /i, is re- 
lated to the area distortion of the lens mapping: 



1/ det Aij = 1/[(1 - ft) 2 - 7 2 



(4) 



where the total shear is 7 = VTi+tI • 

By setting the location of the light ray as the origin of 
coordinates in the lens and source planes, one finds for small 
changes in the position of the light ray: 



(5) 
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The interpretation of this mapping is as follows: if the gra- 
dients of the shear and convergence do not change signif- 
icantly across the image, then the first-order weak lensing 
effect is to produce an additional ellipticity in the shape 
of background sources. While this approach is satisfactory 
for point sources, it breaks down for extended sources. This 
is seen even for the conceptually simplest lens model, the 
point mass or Schwarzschild lens, which produces increas- 
ingly "arcy" images for extended sources at small impact 
parameters to the lens. 

In order to account for the effects of gradients in the 
shear and convergence, a second-order Taylor expansion of 
the gravitational field is required: 



9j + ^D ljk 60 3 Se k . 



(6) 



Following the approach of Bac on et al.l (|2006h , the tensor 
term, Dijk, is associated with two additional lensing distor- 
tion terms 



D 



ijk 



ijk- 



(7) 



By further defining T — T\ -\-iTi and Q — Q\ -\-1G2, we note 
that the components of first (J 7 ) and second (Q) flexion can 



be written: 




T\ = -i(Diii+Di22), 


(8) 


Ti — — — (D211 + -D222), 


(9) 


Gi = —(Dm -3^122), 


(10) 


Q2 — —-(3^211—^222)- 


(11) 



The components of first and second flexion only require a 
subset of the tensor components, Dijk, and do not depend 
on knowledge of D112, D121, ^212 or 1)221- The relevance of 
this is discussed in Section [331 



we can now express the shear, convergence and 
second flexion with their explicit dependence on 
coordinates and surface mass density (see Paper 
details) : 



7i 



72 



T 2 = 



Yj cy 
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first and 
physical 
I for full 



(17) 
(18) 
(19) 
(20) 
(21) 
(22) 
(23) 



3 FLEXION WITH THE RAY-BUNDLE 
METHOD 

In this section, we describe how the ray-bundle method can 
be used as a numerical means of estimating flexion along 
a given line-of-sight. We use the Schwarzschild lens model 
(hereafter, SL), as it provides us with simple analytic so- 
lutions for all of the relevant lensing properties. However, 
since the SL has zero first flexion everywhere except at the 
origin, we emphasise the recovery of second flexion with the 
RBM, whilst demonstrating that results remain consistent 
with zero first flexion. We consider both backwards (lens 
plane to source plane) and forwards (source plane to lens 
plane) ray-bundle propagation, as this separately allows us 
to constrain the appropriate bundle radius to use and to 
quantify the extent of the flexion zone. 



2.3 Circularly symmetric lens models 

For a circularly symmetric mass profile, the deflection angle 
is 



_4G M(|) 



where |£| = Vsi + an d the projected mass, M, is 
rl£l 



M(|e|) = 2^ / £(0£'d£', 
^0 



(12) 



(13) 



for surface density, This enables us to re- write equa- 

tion (pQ) as 



1 



1 M(|C|) 



7rS cr |£| 2 
with critical surface density 

cr AttGDlDls' 
Defining the function 

<2(C):=^S(C)C 2 -M(C) 



(14) 



(15) 



(16) 



3.1 Inverse ray-tracing 

While knowledge of an image location uniquely defines the 
source location, in general, equation (pQ) is not easily in- 
vertible to give all image locations for a given source posi- 
tion. While "brute force" solution methods can be used (e.g. 
Paczyhski 1986), the main alternative is to use inverse ray- 
tracing in its direct form, as was introduced by Kayser, Refs- 
dal & Stabell (1986) and Schneider & Weiss (1986;1987), or 
in its hierachical tree-code form (Wambsganss 1990; 1999). 
Here, light rays are projected backwards from the observer, 
through the lens plane to the source plane, which is rep- 
resented by a two-dimensional grid of source pixels. The 
deflection of each light ray is calculated with equation JTJ). 

While ray-tracing methods are extremely well-suited to 
studying statistical lensing effects (e.g. the creation of mag- 
nification maps for studying probabilities of high magnifica- 
tion events in quasar microlensing) , they are less well-suited 
for studying lensing effects along a given line-of-sight. If it 
were possible to write down an explicit analytic form for 
the null geodesic equation of general relativity for an arbi- 
trary mass (i.e. lens) distribution, the optical scalar equa- 
tions (Sachs 1961; Dyer & Roeder 1974) could be used to 
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Image 



Source 





Figure 1. Shapes in the image and source planes showing non- 
linear lensing effects for extended image/source pairs. For the 
simplest case of a Schwarzschild lens, a circular image is mapped 
to a "teardrop" source (top row), and a circular source is mapped 
to an arc-like image (bottom row). 



measure the changing shape of a (small) bundle of light rays 
as it propagates from the source to the observer. Unfortu- 
nately, specific solutions only exist for a limited number of 
cases, such as light propogation in Swiss Cheese "inhomoge- 
neous" cosmological models (Harper 1991), so an alternative 
approach is required. 



3.2 The ray-bundle method 

The ray-bundle method (RBM) was developed as an alter- 
native to grid-based inverse ray-tracing for studying grav- 
itational lensing along specified sightlines (Fluke, Webster 
& Mortlock 1999). An image shape is denned in the im- 
age plane (or first lens plane in the case of multiple lens 
planes) , comprising a central light ray, which represents the 
"null geodesic" , surrounded by a bundle of iV ray light rays. 
The light rays in this bundle are traced backwards to the 
source plane, producing a distorted source shape. As such, 
RBM provides a numerical analogy to the optical scalar so- 
lution, while retaining properties of the inverse ray-tracing 
approach, such as the use of the thin lens equation for de- 
flection calculations. 

In earlier work (Fluke, Webster & Mortlock 2002), the 
RBM was used to obtain magnification probability distri- 
butions for dark matter-only cosmological models. However, 
the RBM's emphasis on bundle shape means it is ideally 
suited for studying both linear and higher-order gravita- 
tional lensing phenomena (viz. shear and flexion) in detail 
along any given line-of-sight. 

Moreover, the RBM gives a great deal of flexibility in 
choosing the size and shape of the initial bundle (i.e. the 
image) , so that we can constrain the length-scale over which 
the flexion approximation of equation J6} is valid. If the gra- 
dient of the shear and convergence across an image is negli- 
gible, i.e. we are in the regime of equation (0, which is the 
truly weak-field limit of gravitational lensing, then a circular 
source will appear as an elliptical image, and a circular im- 
ag41 would be indicative of an intrinsically elliptical source. 
When tidal gradients do occur across an image, most read- 
ily due to the extended nature of image/source pairs, the 



although unlikely to occur naturally. 



relationship betwen image and source shapes is less obvi- 
ous. For the SL, circular sources are mapped to "arc-like" 
images, and circular images are due to "teardrop" shaped 
sources. This correspondence is shown in Fig. [1] 



3.3 Flexion with the ray-bundle method 

By analogy with equation JS} , small changes in the position 
of a light ray in the image (661,662) and source (S/3i, Sfo) 
planes are related through the following equations, which 
are explicit expansions of equation ©: 



6P1 = 



92 + ^iii^(9? + 



A 11 66 1 +A 1 
+^[D 12 i+Dii2]60i60 2 , 



h + A226O2 + -D21U 



9? + 



1 2 
-D122662 



1 2 
-D222662 



(24) 



(25) 



6f3 2 = A21 

+ ^[D 2 12 + D 2 2l] 661662. 

By choosing the number of rays in each bundle, it is possible 
to obtain a complete solution for the unknown Aij and Dijk 
terms along an arbitrary line-of-sight. Each light ray in the 
bundle experiences the tidal field around the central light 
ray, and so each pair of image and source rays provides a 
unique solution to equations (|24[) and (|25[). Since there are 
five unknowns we wish to solve for, we need only use six 
light rays in total per ray-bundle: a central light ray plus 
^V r ay = 5 rays defining the circumference. Writing these last 



equations in matrix form for each of the five 
bundle pairs (k = 1 ... 5): 

/ 



and 



7ii 

9l2 



0V21 
6622 



7HC 
712(5 



711 
712 



7110021 
?12 0*022 



0(7210(721 
66226622 



0(7i 5 0(7i5 0(7i5 

An 
A12 

\Dm 

~[£>112+£>12l] 

^D122 



/ sen 

66 12 



V 66! 



6621 
6622 



6615661 



\ 



A21 
A22 
\D211 
212 + D221] 
\D222 



825 < 


56256625 


\ 


( SPu \ 




5(3l2 




Spi 3 




SPu. 


/ 


\ sp 15 J 


021 < 


562l6621 


822 < 


56226622 


025 < 


56256625 


\ 


( SP21 \ 




5p22 




8P23 




5p2i 


/ 


\ SP25 J 



(26) 



(27) 



The unknown column vectors of the Aij and Dijk terms 
are the solutions to a set of simultaneous equations. We 
note that there is a degeneracy in the co-effecients of the 
cross-terms, 661662, which means we cannot solve uniquely 
for D112, Di2i, D212 and D221. Instead, these terms can be 
obtained from symmetries of D^k, viz. D^k = Dikj and 
Di2i — D 2 ii- However, as we highlighted earlier, these four 
elements of the Dak tensor do not contribute to first or 
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second flexion [see equations (|8t lll[) ]. and are instead asso- 
ciated with the twist and turn operators identified by Bacon 
& Schafer (2009). 

3.4 Image-to-source plane ray-bundles 

Our first aim is to examine the (physical) scale in the image 
plane over which the RBM can be used to recover flexion. 
Whereas we used angular coordinates for the lens equation 
and flexion terms in sections 2.2 and 3.3, it is now more 
convenient to use a scaling relative to the size of the angular 
(point-mass) Einstein radius: 



0e 



4GM D LS 



(28) 



c 2 DsD L 

in radians, for a lens with mass, M; G and c are the gravi- 
tational constant and speed of light respectively. 

In terms of scaled coordinates, x% — Si/Os and yi = 
A/#e, the two-dimensional lens equation for the SL model 
has solutions: 



Ui — %i 



1 



X -| — \~ X <2 



X% — 



2 



i±Wi + 



2/1+2/2 



2/i + 2/2 / 0, 



(29) 
(30) 



so that the Einstein ring has a radius of 1 unit, and we can 
specify the bundle radius, Ax, as a fraction of the Einstein 
radius. This scaling changes the Taylor expansions: 

r 2 21 

Syi = AuSxi + Ai 2 Sx 2 + — [D 111 Sx 1 + D 12 2Sx 2 \ (31) 
+ — [D 12 i + D 112 ] Sxi5x 2 



E r 2 2~l 

Sy 2 = A 21 Sxi + A 22 Sx 2 + — [D 211 Sx 1 + D 222 Sx 2 \ (32) 



0e 



+ — [D 2 i 2 + D221] Sxi5x 2 



and the column vectors in equations (28) and (29) are now 

/ A ^ \ ( A * \ 

and 



An 
[-D112 +-D121] 

^£>122 



V 



A21 
A 2 2 

°-fD 2 ii 
[-D212 + D221] 

^£>222 



(33) 



once we substitute Syi = S/3i/0E and Sxi = 50i/0E- Note 
that the dimensionless nature of shear, through the Aij 
terms, is maintained with this coordinate change: 



7i 



72 



■ X 1 



—2x\x 2 



(34) 
(35) 

■ Fijk+Gijk (Bacon 
et al. 2006), which reduces to D^k = Gijk for the SL model, 
as all the Fijk = 0. In terms of components, 



as expected. 

Next, we use the relationship D^k 



D 



ijl 



D 



ij2 



— Giji — - 

— Gij2 — + 



Gi 

G2 



-Gi 

Gi 
G2 



(36) 
(37) 



and we use equations (|22|) and (|23|) , to calculate the analytic 
solutions for the two components of second flexion for the 
SL model: 



1 4xi [x\ - 3x|) 
" fe (xl + x 2 2 f ' 

1 4x 2 (3xj - x|) 
~ {xi+xlf 
Combining these last two equations implies 

\G\= 1 ' 4 



2\3/2' 



(38) 
(39) 

(40) 



Fortuitously, the O^ 1 term in equations (|38|) and (|39p ef- 
fectively cancels the #e coefficient of the D^k terms in the 
column vectors, equation ((33]), so we can obtain the scale in- 
dependent second flexion with the RBM without requiring 
a specific value of #e- 

Each ray-bundle is characterised by the radius, Ax, two- 
dimensional location of the central ray, Xi, and the number 
of bundle-rays, iVray As shown in Section 3.3, we choose 
^V ray = 5 in order to solve for the unknown Aij and Dijk 
terms. For convenience in plotting and making comparisons 
with analytic results, we use polar coordinates, (r, <fi), where 
r — \J x\ + x\ is the radial impact parameter, and cj) — 
tan _1 (x2/xi) is the polar angle. The analytic and RBM- 
recovered solutions are sampled on a grid with dimension 
N r x in polar space, which means that the sampling of 
solution space is sparser with increasing r. 

For each Ax, we have to choose an appropriate range 
of r values for sampling. The minimum impact parameter 
is r m in = 1 + Ax(l + e), where e — »• is a small numerical 
offset (we use e — 0.01). The offset avoids light rays in the 
image bundle overlapping the Einstein radius, thus ensur- 
ing that the 5x5 square matrix of 5x terms is invertible. 
While the lens equation describes a mapping from the image 
plane backwards to the source plane, not all RBM images 
are permissable. For example, it is not possible to have a 
regular polygonal image straddling the Einstein ring that 
corresponds to a single source shape. However, the RBM 
does work for images either completely outside the Einstein 
radius, or completely within it^For more complex lens mod- 
els, where there are non-degenerate caustics, the solution is 
less obvious. 

The upper limit, r ma x, for the impact parameter is cho- 
sen to be Max (10, lOAx), with the latter limit only relevant 
for bundle radii greater than the Einstein radius. In prac- 
tice, the majority of bundle radii we use were Ax < 0.1 (see 
Fig. 6). Our choice of r max = 10 was based on numerical 
tests, which included varying the number of radial samples, 
N r . Beyond this impact parameter, differences between the 
analytic and RBM-recovered flexion values were consistent 
with numerical (i.e. precision) limits. 

For each bundle ray, k = 1 . . . A^ ray , we determine 
(5xik,Sx 2 k) relative to the central ray, and either apply 
equations (|3T|) and ((32}, assumed to be the exact analytic 
solution, or utilise the RBM approach to deflect light rays 
with the gravitational lens equation, to obtain (Syi, Sy 2 ). We 

3 Images inside the Einstein radius are in the strong lensing 
regime, and hence outside the flexion zone that we define in Sec- 
tion (331 
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Figure 2. Qi(r, </>)#e for an image bundle radius Ax = 0.01, 

r = \jx\ + x\ is the radial impact parameter in the range 1 + 

Ax(l + 0.01) < r < 5.0, (j> = tcni- 1 (x 2 /x 1 ) and x = 0/0 E . The 
surface is sampled on an 100 x 100 polar grid. 



then build the matrices of simultaneous equations for each 
ray-bundle and use Gaussian elimination with back substi- 
tution to solve for the Aij and Dijk terms, and hence G% and 

To test the matrix inversion code, we used the SL so- 
lutions of equations (|34H35[) and (|38H39[) . as the analytic 
solutions for each (r, 0) sample. Within the limits of nu- 
merical accuracy afforded by our implementation, we find 
that we correctly recover the input Ti = and Qi, inde- 
pendent of N r , and Ax. This demonstrates the accuracy 
and utility of our simultaneous equation-solving code, under 
the proviso that there is no degeneracy in the mapping of 
image ray-bundles to source ray-bundles. This assumption 
is appropriate in the "almost weak" regime where flexion 
acts, and by ensuring that no part of the source crosses the 
caustic point for the Schwarzschild lens. 

Next, we use the scaled lens equation, equation ([29]) . to 
deflect the individual light rays in each bundle, so that we 
have pairs of image and source bundle shapes, and build the 
matrices for each bundle. Inverting the matrices, we obtain 
RBM-estimates for second flexion and use recovered first 
flexion results as a consistency check on numerical effects 
(see below). 

In Fig. [2j we plot an indicative second flexion surface, 
Qi(r, 0)0e; second flexion values must be scaled by the an- 
gular Einstein radius (in radians) for comparison with a 
particular system. The spin-3 oscilliatory nature of second 
flexion is visible, with the amplitude increasing towards the 
origin. The surface is plotted over a limited range in im- 
pact parameter (1.0101 < r < 5.0), for an image bun- 
dle with radius Ax = 0.01. The surface is sampled over 
& N r x N(j> = 100 x 100 polar grid. 

We compare our RBM-based estimates of second flexion 
with the analytic values using the mean-square error, M s : 




Figure 3. Peak signal-to-noise ratio (Ps) versus image size, Ax, 
for RBM-recovery of Q\ . 



N X N< 



N i N 2 

^^|G(p,4)-6(p,</)| 2 . 



(41) 



p=l q =l 



Here, G(p, q) is the grid-sampled analytic surface (one of Qi 
or |(J|), and G(p,q) is the grid-sampled, RBM estimate. We 
define a peak signal-to-noise ratio as: 



Ps = 201og 10 



Max(G) 



(42) 



which provides a quantitative value for the equality of sur- 
faces: Ps — »• oo if the surfaces are identical, and Max(Gr) is 
the maximum value of the surface G(p, q). We find minimal 
dependence in the calculated Ps on the gridded surface res- 
olution, with a variation APs — ±0.4 when Ps > 20 for 
grids with dimensions 50 < N r , < 250. 

We plot Ps versus Ax for RBM-recovery of Qi in Fig. 
[3] Results for Q2 and \Q\ are comparable with Q\. For the 
components of second flexion, the peak signal-to-noise con- 
tinues to grow as image size is reduced, and we approach 
the theoretical point-source/point-image case. This result 
demonstrates that the RBM can indeed be used to recover 
second flexion for the SL model, and gives us confidence that 
this technique can also be applied for other, more complex 
lens models in a backwards ray-tracing mode (i.e. when the 
lens equation cannot be inverted to give image locations as 
a function of a source location). 

Since the SL model has zero first flexion, we are not able 
to demonstrate that the RBM is actually sensitive to the 
centroid shifts from a lens model with non-zero first flexion 
components. What we can show, however, is that RBM is 
consistent with recovering zero first flexion for the SL model. 
Using an approach similar to that outlined above for Qi and 

we determine estimates for T\ from the matrix solutions 
as shown in Table [1] The maximum (absolute) values of the 
recovered T\ values are based on bundle radii in the range 
10 -4 < Ax < 1 and polar grid resolution 50 < N r ,Ncp < 
500. Results for Ti and \T\ are comparable. 
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Table 1. Limits on recovery of first flexion, expected to be T\ = 
0, for the Schwarzschild lens model using the RBM. Ax is the 
image bundle radius in units of the Einstein radius. 



Ax 


Max(|^^ E |) 


0.0001 


1.1 x 10~ 7 


0.001 


2.2 x 10~ 9 


0.01 


1.9 x 10~ 6 


0.1 


1.1 x 10- 3 


1.0 


3.1 x 10~ 2 



c c 




F F 



Figure 4. Pairs of sources and images represented as ray-bundles 
with N ray = 5 light rays around a central ray for circular images 
(thin lines) and circular source (thick lines). The dotted line rep- 
resents the Einstein radius for a Schwarzschild lens. 

A slight increase in the maximum recovered first flex- 
ion is seen at the lowest bundle radius: from ~ 10 -7 for 
Ax = 10~ 4 to - 10~ 9 for Ax = 10~ 3 , suggesting that we 
have reached a numerical limit (in our implementation) . We 
propose that Ax = 10 -3 is thus an appropriate lower bun- 
dle radius to maximise accuracy in using the RBM for an 
arbitrary lens configuration for recovering flexion. 

Due to our choice of A^ ray = 5, there is a slight depen- 
dence on the orientation of the bundle with respect to the 
lens. Suppose we create bundles that are evenly spaced in 
angle, 0, around the lens, so that the bundle centres are: 

I c (r, <j>) = [xi,x 2 ] = [r cos(0),r sin(0)] , (43) 

and the individual bundle rays are located at two- 
dimensional coordinates 

I fc (r, 0) = [xi + Axcos(^fc + ^o),x 2 + Axsin^fc + ^ )] , (44) 

with ipk = 2k7v/N ray for k = 1 . . . iV ra y, and ipo is a phase 
term that is fixed per bundle. 

If we set ^0 = for all bundles, we get the situation 
shown in the left-hand panel of Fig. [4) where the image (reg- 
ular pentagons) and source bundles (irregular pentagons) 
are shown with respect to the Einstein radius (dotted line). 
While the bundles retain the same orientation in the im- 
age plane, independent of 0, the radial separation between 
the lens and the closest point of the bundle varies with <\>. 
This means that the tidal field experienced by each bundle 
depends weakly on 0. For the case of image bundles that 
appear more circular (right-hand panel, iV ra y = 100) this 
effect is much less pronounced, and we have the expected 
rotational symmetry. 




Figure 6. Maximum (upper, thick line) and minimum (lower, 
thin line) percentage error, r% in RBM-recovered Q(r) as a func- 
tion of image bundle radius, Ax. 

The consequence of this is demonstrated in Fig.[5j where 
we plot \G(r, 4>)\ on a 100 x 100 polar grid for Ax = 0.2 and 
^o = (top left), Ax = 0.2 and uniform randomly selected 
ipo (top right), Ax = 0.01 and uniform randomly selected 
(bottom left) and analytic \G(r, <p)\ (bottom right). A sinu- 
soidal variation with (j) is apparent when ip = for fixed r, 
and a scatter is introduced when ipo is chosen uniformly ran- 
domly between [0, 2ir). The magnitude of the (j) dependence 
is reduced by choosing a smaller Ax. 

To quantify this effect we determine the minimum im- 
pact parameter, r T , at which ^G(rv) — G(r r )j /G(r) < r% 

for a given image size, and G(r T ) — Max [G(r, (/>)]. In Fig. [6] 
we plot the maximum (thick line) and minimum (thin line) 
percentage deviation between |<?(r)| and |<?(r)| as a func- 
tion of the image bundle radius. Samples were evaluated on 
a 250 x 250 polar grid. The horizontal line is a nominal error 
of 1%, indicating that a sufficiently high level of accuracy 
occurs for Ax = 0.01 (vertical line). Moreover, we find that 
for Ax < 0.01, the error in G% is less than 1% for all impact 
parameters r > 1.01, which suggests that we can use the 
RBM approach for impact parameters extending almost to 
the Einstein radius. From this result, we infer that any vari- 
ation we see in the RBM-recovered G% for an arbitrary lens 
configuration is evidence of a ^-dependent second flexion, 
which differs from the predictions of a circularly-symmetric 
lens models. 



3.5 Source-to-image plane ray-bundles 

For the SL we have the fortunate situation that there are 
analytic solutions for both 0^/3 (backwards) and /3 
(forwards) versions of the lens equation. This gives us more 
flexibility in denning the image shape in order to test the 
range over which the flexion formalism is appropriate. In the 
previous section, we started with a regular image, which is 
deflected to a teardrop source. Here, we use a regular polygo- 
nal source, and obtain the corresponding image shape. Prop- 
agating this modified (i.e. flexed) image backwards as a ray- 
bundle, we hope to recover the input source shape. 

We note that an image can cover a small portion of the 
image plane, and hence the tidal fields can be weak over 
the extent of the image. For a strongly flexed ("banana- 
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Figure 5. Sample \G(r, 0)|#e surfaces obtained with the ray-bundle method, demonstrating the effect of the phase angle, ipo, for bundles 
in the image plane. Each surface is calculated on a 100 x 100 polar grid, with impact parameters in the range 1 + Ax(l + 0.01) < r < 5. 
(Top left) Recovered with Ax = 0.2, bundle phase i/jq = 0. (Top right) Recovered with Ax = 0.2, random i/jq. (Bottom left) Recovered 
with Ax = 0.01, random ipo. (Bottom right) Expected analytic \Q(r : </>)|#e- 



shaped") image, more of the image plane is covered by the 
image, and hence tidal effects are more prevalent. For an 
arbitrary lens model or configuration, where we do not have 
analytic solutions for f3 —> 0, we must work in the "regular 
polygonal image maps to unknown source shape" regime, so 
it is critical to understand the limits of our technique. 

Fig. [7] shows six sample configurations. In each case, a 
circular source (A) is mapped to a flexed image (B) using 
the inverse solution, equation ([30]) . of the lens equation for 
the SL. Next, equation © maps the differential image ray 
coordinates, S0i, back to the source plane to produce a new 
source profile (C). Finally, equation J5} is used to create 
an elliptical, shear-only image shape (D) from the circular 



source. The major and minor axis lengths of image D are 
comparable to the equivalent (distorted) axes for the flexed 
image (C). Source C and image D are presented as mir- 
ror images of their true positions for clarity in the figure. 
We have used N ray = 100 for the bundles in this figure in 
order to show the true shapes, hence we can neglect the <fi- 
dependence (see Section 1X4)) . While the Scwharzschild lens 
model produces two images for each source position, we only 
consider the more highly magnified image outside of the Ein- 
stein radius. Significant variation between the source shapes, 
A and C, is an indication that we are in a regime where the 
second order Taylor series expression for flexion is not valid 
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Figure 7. Sample plots showing the limitations of the flexion formalism for extended sources using the Schwarzschild lens model, with 
respect to the Einstein radius (dotted circles). In each case, a circular source (A) is mapped as a ray-bundle to a flexed image (B) 
using the lens equation (i.e. reverse RBM). This flexed image is mapped back to the source plane (C) using equation (0. The elliptical, 
shear-only version of the image (D), obtained by applying equation ((5} to the circular source, is shown for comparison. Note that source 
C and image D are presented as mirror images of their actual positions for clarity in this figure. Analytic expressions for flexion of 
extended sources are valid in cases where source shapes A and C are comparable; shear-only weak lensing analysis is satisfactory when 
image shapes B and D are comparable. The level of discrepancy between true and recovered source shapes increases for larger image 
radii and small impact parameters. Each bundle comprises 100 rays. Source bundle radii are: Ay = 0.5 (left column), 0.2 (top right and 
middle right), and 1.0 (bottom right). X and S values are noted for each scenario, as defined in equations (|45|) and (|46|) . 



for extended sources - this is the true strong lensing regime 
where images are arcs rather than arclets. 

We examine the level of agreement between the Taylor 
expansion, equation J6}, and the forwards RBM solution by 
comparing the relative locations of each image or source ray 
in the bundle. We introduce two quantities: 

X = log 10 aj = log 10 I j^— l 1 ^ - l' D ,n\ 2 J , (45) 

where the two-dimensional vectors, 1b, n and Ii> jTl , are the 
iV r ay light rays in the actual (B) and elliptical (D) images, 
and similarly 

/ 1 ^ray \ 

S = \og 1Q * 2 s = log 10 I — \ Sa ^ ~ S 'c,n\ 2 J , (46) 

where the two-dimensional vectors, Sa,u and Scr n , are the 
light rays in the initial (A) and recovered (C) source bundles. 



We calculate X and S for source bundles with radii Ay = 
0.01,0.02,0.05,0.1,0.2,0.5, 1 and 2. In Fig. [3 we plot X and 
S as functions of the image bundle impact parameter, x c . We 
select X > — 4 and S > —4 as indicative that the distortion 
is significant; this was confirmed by eye using plots similar 
to Fig. [7] For higher values of X and <S, there were clear 
differences between sources A and C, and images B and D. 
We find 

x c > 1 + 2.2 Ay, (47) 

for comparison of source bundle shapes A and C (<S > —4) 
and 

x c > 1 + 3.6Ay, (48) 

for comparison of image bundle shapes B and D (X > —4). 

Next, we consider a more quantitative approach, based 
on comparing bundle ellipticities, in a manner comparable 
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-16 L , , , , , J -16 L , , , , , ^ 

-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 

logiofe- 1 ) l°g 10 (* c -l) 

Figure 8. X (left) and S (right), as denned in equations (|45|) and (|46|) , plotted as a function of log 10 (x c — 1). Individual lines are 
for different source radii (top to bottom) Ay = 2.0,1.0,0.5,0.2,0.1,0.05,0.02 and 0.01. The horizontal dashed line in each panel is at 
I = S = —4, taken as the limit above which shape distortion is significant. The dots represent the empirical fits (left) x^ > 1 + 3.6Ay 
for X and (right) x\ Q > 1 + 2.2 Ay for S. The noise in the lowest lines is due to the use of single precision floating point numbers. 



to the standard analysis for examining weak lensing-induced 
shear. Defining quadrupole terms 



I 



i^X'i X j ) i^X j 2/j)d X 



relative to the bundle centroid 

iVray 

%i — ^ ^ %i,m 
m=l 



(49) 



(50) 



we consider a complex ellipticity (e.g. Schneider 2005) of the 
form: 

= Qii - Q22 + 2iQ 12 
Qn + Q22 

which has norm: 



X 



1x1 = 



y/(Qn -Q 22 ) 2 + 4Q 



Qn + Q22 

We define error terms: 
Eac = IxU ~ \x\c aiul 



Ebd 



IxU 

\x\b - \x\d 
\x\b 



(51) 
(52) 

(53) 
(54) 



for comparison between ellipticities, and determine the im- 
age plane impact parameter, x c , at which a source with ra- 
dius, Ay, first exceeds E — 1%, 5% and 10%. 

Since a circular source has |x| = 0, we use elliptical 
sources with axis ratios b/a — 0.8, 0.9 and 0.99. We consider 
the two cases where the semi-major axis is aligned tangen- 
tially or radially to the Einstein radius, computing these 
limits and also the average (based on original data values) 
of these two orientations. 

Plotting results as log 10 (x c — 1) versus log 10 Ay, see 




-1.5 -1.0 -0.5 

log 10 (Ay) 



Figure 9. Minimum impact parameter, x c , at which the error, 
EaCi is first below 5% as a function of the source radius, Ay, 
for elliptical source shapes with axis ratios b/a = 0.8 (red, dot- 
dashed line), 0.9 (blue, dashed line) and 0.99 (black, solid line). 
Lines are the least-squares fits in log-log space to the functional 
form x c = 1 + eAy n . 



Fig. [9l we see a relationship that is highly suggestive of a 
functional form: 



1 + eAy n 



(55) 



We perform a least-squares fit in log-log space to obtain the 
parameters e and n. Results of these fits are presented in 
Table [2] - in all cases, the calculated Pearson coefficient is 
r > 0.994, indicative that equation (|55p is an appropriate 
functional form. There is variation in the fitted parameters 
based on the chosen source axis ratio. This is not suprising, 
as the tidal gravitational field across the resultant image 
depends on the relative separations and orientation of indi- 
vidual image rays from the lens (c.f. with discussion on ori- 
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Ay 

Figure 10. The boundaries of the flexion zone (blue line) and 
the shear zone (red line), as defined by equations (|56|) and (1570 
respectively. The white region is the true strong lensing regime, 
while the light grey-shaded region represents the preferred zone 
for both weak lensing shear and flexion analysis. The dark grey- 
shaded region, Ay > 0.5 and hence source sizes comparable to 
the Einstein radius, means the stronger tidal fields "flex" rather 
than "stretch" images. 



Table 2. Least-squares fitting parameters for the functional form 
x c = 1 + eAy n based on relative errors in source A and C elliptic- 
ities. Fits are made for sources with axis ratios b/a = 0.8, 0.9 and 
0.99 with the semi-major axis aligned either tangentially or radi- 
ally to the Einstein radius. Fits were performed in log-log space, 
but averaging is performed with original data values. In all cases, 
the Pearson coefficient is r > 0.994. The inferred e and n values 
based on the limit *S > —4 are shown for comparison. 



Tangential Radial Average 



b/a 


Eac 


e 


n 


e 


n 


e 


n 




1% 


3.58 


0.65 


2.14 


0.56 


2.86 


0.61 


0.8 


5% 


2.27 


0.67 


1.60 


0.63 


1.94 


0.65 




10% 


1.87 


0.68 


1.35 


0.65 


1.61 


0.67 




1% 


4.05 


0.61 


2.95 


0.55 


3.50 


0.59 


0.9 


5% 


2.69 


0.65 


2.16 


0.62 


2.43 


0.63 




10% 


2.25 


0.66 


1.82 


0.64 


2.04 


0.65 




1% 


5.67 


0.50 


5.08 


0.48 


5.38 


0.49 


0.99 


5% 


4.32 


0.56 


4.08 


0.59 


4.20 


0.56 




10% 


3.77 


0.58 


3.50 


0.58 


3.63 


0.58 


Average 


5% 










2.70 


0.62 




S 










2.2 


1.0 



Table 3. Least-squares fitting parameters for the functional form 
x = 1 + eAy n based on relative errors in image B and D ellip- 
ticities. In all cases, the Pearson coefficient is r > 0.995. The 
inferred e and n values based on the limit X > —4 are shown for 
comparison. 



E BD 


e 


n 


1% 


3.87 


1.15 


5% 


1.76 


1.17 


10% 


1.20 


1.15 


X 


3.6 


1.0 



second-order Taylor series approximation given by equation 
((6]) is not sufficiently accurate when applied to an extended 
source, and we are in the true strong lensing regime. The ra- 
dius of the flexion zone boundary increases for larger source 
sizes, Ay, relative to the Einstein radius, which is expected 
as there will be greater tidal field variations across an im- 
age/source bundle. 

We perform a similar analysis for the relative elliptic- 
ity error between images B (RBM) and D (shear-only), al- 
though we now revert to using a circular source only. Pa- 
rameters are presented in Table [3] We find that for 

x c > 1 + 3.87A2/ 1 - 15 , (57) 

the lensed image shape and a shear-only intereptation are 
essentially the same (error < 1% in ellipticity) implying 
the first-order Taylor expansion is adequate. We refer to 
the boundary defined by this expression as the start of the 
"shear zone", where image shapes are essentially indistin- 
guishable from ellipses, and a traditional weak-lensing (i.e. 
shear and convergence only) analysis is satisfactory. How- 
ever, our result does not preclude use of the flexion formal- 
ism at larger impact parameters, as a potentially measurable 
non-zero flexion remains. 

There is a crossover between the boundaries when 
Ay ^0.5, indicated by the dashed line in Fig. [TUJ In some 
sense, the light grey-shaded region is the preferred region for 
weak lensing shear and flexion analysis, as it corresponds to 
both shear-only ellipticity errors < 1% and flexion-recovered 
source ellipticity errors < 5%. The dark grey-shaded region, 
Ay ^0.5 and hence source sizes comparable to the Einstein 
radius, means the stronger tidal fields "flex" rather than 
"stretch" images. The effect of this is demonstrated more 
clearly in Fig. 12, when we consider cosmologically-realistic 
scenarios. 



entation of image bundles in section 3.4). As a best estimate, 
we average over the three chosen axis ratios for Eac = 5%, 
to obtain the second-to- last row in Table [2] 

We propose the following interpretation: if we see an 
image that looks like B, and we use the flexion formalism 
to determine what the source would look like, we would be 
wrong (error of > 5% in source ellipticity) unless: 

x c > 1 + 2.7A2/ - 62 . (56) 

We refer to the boundary defined by this expression as the 
start of the "flexion zone" (blue line and both grey-shaded 
regions m Fig. [10]). For image impact parameters closer 
to the lens than this limit (white region in Fig. I10p , the 



4 APPLICATION 

We approach our application of the results from the preced- 
ing section with the understanding that the Schwarzschild 
lens is not an ideal description of the extended mass distri- 
bution of, for example, a galaxy cluster lens. Analytic pre- 
dictions for first and second flexion do depend on the density 
profile of the lens model. Indeed, from Fig. 2 in Paper I, it 
can be seen that the Schwarzschild lens has zero k and 
while these values are non-zero for the extended mass pro- 
files (ie. SIS, NFW and Sersic profiles). However, BirkhofT's 
theorem allows us to consider all truncated lens models to 
be Schwarzschild- like outside of the truncation radius. 
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Figure 11. The RBM-based constraint on bundle sizes for re- 
covery of analytic flexion results applied to extended sources. For 
lens masses M/M = 10 13 , 10 14 and 10 15 , we plot the quantity 
f° r fixed Zd = 0.1 (top) and Zd = 0.5 (bottom), as a func- 
tion of z s , with Zd < z s < 2. Results are for Ar] = 10 kpc. The 
lower dashed line in each panel indicates 1% errors in recovered 
flexion, and the upper dashed line is for a 10% error. 



smaller than the radius of a fiducial Milky Way-sized galaxy 
(Arj = 15 kpc at z s =0), noting that the curves scale lin- 
early with bundle radius. We use the concordance cosmology, 
with total matter density, ^m,o = 0.3, dark energy density, 
^a,o = 0.7 and Hubble parameter, Ho — lOO/i km s _1 with 
h = 0.7. 

The desired criteria (lower dashed line) at 0.01 for a 1% 
flexion error, is only met for low lens redshifts and high mass 
lenses, typical of galaxy clusters, and for small source galaxy 
radii. A more conservative limit at 0.1 (upper dashed line), 
results in a 10% error in recovered flexion values for extended 
sources - see Fig. [6] While this may appear somewhat dis- 
couraging for flexion programs, we note that the required 
size criteria could be reached by considering isophotes of 
an image, corresponding to smaller source sizes. We do not 
discuss this further in the present work. 

In Fig. 1121 we plot the location of the inner boundaries 
of the flexion zone (dashed lines) and the shear zone (solid 
lines), for several typical scenarios: lens masses 
(top panel), 10 14 M o (middle panel) and 1O 15 M (bottom 
panel). In each panel (from top to bottom), the lens red- 
shifts are Zd — 0.1 (black), 0.2 (red), 0.5 (green) and 1.0 
(blue), the source radius is A77 = 10 kpc, and Zd < z s < 2. 
Note the crossover between these regions that occurs for 
M = 1O 13 M , and in most cases presented here, the shear 
zone actually starts closer to the lens than the flexion zone. 
Fig. [13] shows the effect of changing the source radius, 
Zd = 0.2, with A 77 = 20 kpc (dashed lines) and 5 kpc (solid 
lines). Lens masses used were 1O 13 M (black), 1O 14 M (red) 
and 1O 15 M (blue). 



For comparisons with observations, we need to convert 
our Einstein-radius scaled results from Section [3] back to 
angular units on the sky. With 0i as the angular position 
in the lens plane and ft the angular position in the source 
plane, we have in the small angle limit (which is appropriate 
for sources and lenses at cosmological distances): 



0i = XiO^ and ft = yiO^ 



(58) 



where 0e, the angular (point mass) Einstein radius, was de- 
fined in equation (|28p. 

In the absence of lensing, a source at ft would be ob- 
served at 0i , since from equation © , ft = 0i , and the angu- 
lar extent of the bundle satisfies A/3 = AO. Substituting for 
A/3 = Ar]/Ds, equations (|56p and (|57p now become: 



1 + e 



Ar) 



0e 



(59) 

(60) 
(61) 



The RBM-limit on bundle sizes in angular units implies 
that the analytic flexion terms are most accurate for im- 
ages with radii AO < 0.01#e- We present results in Fig. 1111 
For lens masses M/M = 10 13 , 10 14 and 10 15 , we plot the 
quantity e ^ for fixed Zd = 0.1 (top) and Zd = 0.5 (bot- 
tom), as a function of z s , with Zd < z s < 2, and assume 
that AO ~ A/3. We choose Ar] = 10 kpc, which is slightly 



5 CONCLUDING REMARKS 

While weak lensing with shear is now well-established both 
theoretically and observationally in the galaxy-galaxy lens- 
ing and cosmic shear cases, weak lensing via flexion is still 
in its infancy. 

In Paper I, we considered analytic solutions for 
circularly- symmetric lens models. In this paper, we have 
demonstrated how the ray-bundle method can be used 
to recover the analytic second flexion results for the 
Schwarzschild lens model to high accuarcy, and is consistent 
with the recovery of zero first flexion. Indeed, we recover the 
Schwarzschild model second flexion solutions with errors no 
worse than 1% if bundle radii are AO < 0.01#e- In such cir- 
cumstances, the second-order Taylor series expansion used 
by Bacon et al. (2006) is appropriate for extended sources. 

Furthermore, we have identified the existence of a flex- 
ion zone in the image plane, which can be considered an 
optimal region for applying the analytic flexion formalism 
to extended sources. 

The ray-bundle method now provides us with a valuable 
numerical tool for studying flexion due to complex lens dis- 
tributions, such as asymmetric lens models or cosmological 
structures, where no such analytic solutions exist. 
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Figure 12. The location of the inner boundaries of the flexion 
(dashed lines) and shear zones (solid lines) for several typical 
scenarios. Lens masses are M = 1O 13 M0 (top row), 1O 14 M0 
(middle row), 1O 15 M0 (bottom row). Lens redshifts in each panel 
are (from top to bottom) Zd = 0.1 (black), 0.2 (red), 0.5 (green) 
and 1.0 (blue). Source radius is Ar] = 10 kpc. In most cases 
presented here, the shear zone commences closer to the lens than 
the flexion zone. 
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lines are the inner boundaries of the flexion zone for A77 = 20 kpc 
(dashed lines) and 5 kpc (solid lines). Lens masses are 1O 13 M0 
(bottom, black), 10 14 M o (middle, red) and 1O 15 M (top, blue). 
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